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Abstract 

Nonequilibrium electron transport through the Kondo impurity is investigated numerically for 
the system with twenty conduction-electron levels. The electron current under finite voltage 
drop is calculated in terms of the 'conductance viewed as transmission' picture proposed by 
Landauer. Here, we take into account the full transmission processes of both the many-body 
correlation and the hybridization amplitude up to infinite order. Our results demonstrate, 
for instance, how the exact solution of the differential conductance by Schiller and Hersh- 
field obtained at the Toulouse point becomes deformed by more realistic interactions. The 
differential-conductance-peak height is suppressed below e'^/h with the width hardly changed 
through reducing the Kondo coupling from the Toulouse point, whereas it is kept unchanged 
by further increase of the coupling. We calculated the nonequilibrium local Green function as 
well. This clarifies the spectral property of the Kondo impurity driven far from equilibrium. 

PACS codes/keywords: 72.15.Qm (Scattering mechanisms and Kondo effect), 82.20.Mj 
(Nonequilibrium kinetics), 75.40.Mg (Numerical simulation studies) 



1 Introduction 

Recently, in fine semiconductor devices, the Kondo effect has been observed at very low tem- 
peratures |]l], H, 1^. The indications of the Kondo effect are drawn mainly from the electron- 
transport measurements under various controlled system parameters such as the temperature, 
the gate voltage, the hybridization amplitude, the magnetic field and the bias voltage. In par- 
ticular, the differential conductance dI{V)/dV measured over a wide range of the bias voltage 
V shows a distinctive indication of the Kondo effect: The differential conductance shows a peak 
around the zero bias {V = 0) with the width comparable to the Kondo temperature. The peak 
height grows as the temperature is decreased. These features had been observed in various 
other experimental realizations such as tunnel metallic junctions and metallic point contacts 
@7 H IIIj although in these systems, the system parameters are not so freely tunable as in the 
semiconductor devises. Nevertheless, the series of experiments has been stimulating theoretical 
investigations so far g |, 0, ^ |l|, 0, [||l|, |^, [||l| |l|, |2|||, ||, || H, §, ~ 



The essence of the differential-conductance feature is attributed mainly to the shape of the local 



density of states at the impurity. Note that the density of states contributes to the transmis- 
sion probabihty in the intermediate stage of the transmission processes. It is well known that 
the density of states at the Kondo impurity is of a peak structure with the width comparable 
to the Kondo temperature. This picture yields comprehensive qualitative understanding of 
experimental observations. 

Apart from the practical interest such as to explain the essential features of experimental 
observations, the subject may cast a rather fundamental problem where the many-body cor- 
relation and the nonequilibrium-driving force coexist, and are both important. Without the 
biquadratic correlation, we can calculate the nonequilibrium transport exactly with the 
Keldysh approach On the other hand, the transport coefficient for correlated system is 

given by the Kubo formula, although the application is restricted in the vicinity of equilibrium. 
There have not yet been found any frameworks to command the nonequilibrium transport and 
the biquadratic correlation simultaneously. In fact, in the above mentioned theories, a number 
of perturbative approaches — either the voltage drop or the coupling to the leads are supposed 
to be infinitesimal — are employed. Recently, Schiller and Hershfield, however, succeeded in 
treating the situation rigorously ||2l|, To the best of our knowledge, their solution is the 
first exact solution of a nonequilibrium correlated system. They found that the nonequilibrium 
Kondo model becomes integrable at the so-called Toulouse point where the transverse 
and longitudinal Kondo couplings to both leads are tuned very carefully so as to satisfy the 
Emery-Kivelson condition Their result at zero temperature shows that the differential 

conductance is of a pure Lorentzian shape. This result might be rather disappointing, because 
the pure Lorentzian form had been known to be realized in simple quadratic system |p9| . This 
rather uncharacteristic feature may be due to the fact that the system at the integrable point is 
nothing but quadratic. In fact, in a subsequent paper they considered some perturbations 
to the Toulouse point, and tried to show how the integrable-point result is deformed by these 
perturbations. 

Here, we formulate the nonequilibrium transport phenomenon with many-body correlation 
in a manner that is suitable for numerical simulation. Based on the formulation, we investigate 
the transport through the Kondo impurity. The transport coefficient is evaluated in terms of 



the 'conductance viewed as the transmission' picture proposed by Landauer |3^, 34 . Here, we 
take into account the full transmission processes of both the many-body correlation and the 
hybridization amplitude up to infinite order. The rest of this paper is organized as follows. In 
the next section, we explain the model Hamiltonian, the so-called resonance level model, which 
is equivalent to the Kondo model [^5[ |3^. The equivalence is explained, and the advantages 
to investigate the former are explicated. In Section ^, the formulation of the nonequilibrium 
transport phenomenon, mentioned above briefly, is presented in detail. We also propose the 
algorithm to calculate the nonequilibrium Green function. In Section ^, we give the results of 
numerical simulation. In the last section, we give summary of this paper. 



2 Resonance level model 

In this section, we explain the so-called resonance level model. We simulated this model rather 
than the Kondo model. Although both models are equivalent |3^, it is much more advan- 
tageous to treat the former model rather than the latter. The Hamiltonian of the resonance 
level model is given by, 

^ = E + ^) ^i^fc + E ('^F^ " ^) ^I^fe + + S-c^) + US^"^^^. (1) 

The operator l]^ (^fe) creates a conduction electron of wave number fc in the left (right) lead. 
The wave number fc is distributed uniformly over the range fc = — fcc ^ /cc (^c — tt); here, the 
wave number is indexed in such a way as v^ki — 2i/ {N — 1) — 2/ {N — 1) — 1 = 1, 2, • • • , N). 
That is, there are N conduction-electron levels in each lead. The parameter vpkc{= ujc) is 
regarded as the unit of energy throughout this paper; namely, vpkc = 1. The operator 



creates an electron at the impurity position; that is, 




(2) 



The operators S*^ and are the conventional S = 1/2 spin operators. Hence, the parameters 
Vh and U denote the hybridization amplitude and the density-density correlation strength to 
the conduction-electron bands, respectively. Because we study the nonequilibrium situation 
where both leads are kept half-filled, the parameter eV give the chemical potential difference 
between two leads. The way how we realize such nonequilibrium situation in computer will be 
shown in the next section. 

It is shown that the above model (0) is equivalent to the anisotropic Kondo Model by means 
of the bosonization technique |3^ . The anisotropic Kondo model is given by, 

-Hk = E ("F^+^)4-^'«-+E 

ka ^ ^ ka ^ ' 

+j|l^^flfl_fifi + J^(5+cjc| + ^-cjq), (3) 

where the index a denotes the spin index; a =\ or |. According to the bosonization analysis, 
the following mapping relations should hold p^ ; 

PkJ±^^C, ^_l/2 

y = — ^ — (pfc^c) ' , (4) 

and, 

t/= (l-\/2|l-pfcJ||/2|)/2pfc, (5) 

with the density of states pk = (27rwF)~-'^. The mapping relations tell the foUowings. First, 
at J|| = 2(1 — l/y/2)/pk, there is a very special point at which the parameter U vanishes. 
(Note, however, that quantitative reliability of these mapping relations is skeptical. Precise 
connection between these models is somewhat lost in the bosonization procedure. Therefore, for 
instance, we cannot tell definitely at which point (Vh-U) the isotropic Kondo point is realized.) 
At this point, the resonance revel model ([|) becomes quadratic, and this point, the so-called 
Toulouse point , corresponds to the strong-coupling fixed point of the Kondo model. Second, 
the reduction of the Kondo coupling from the Toulouse point corresponds to the attractive 
biquadratic correlation U < 0, whereas the furthermore increase of the Kondo coupling is 
viewed as the repulsive correlation U > 0. Roughly speaking, the previous rigorous analysis 
is concerned in the Toulouse point j2^, In the analysis, however, the authors used the 
model with two channel density-density-correlation term, and arrived at the Toulouse point of 
the two channel Kondo model. Because our simulation is rather incapable of the two-channel 
Kondo physics, we have chosen the model with the conventional Kondo coupling. We investigate 
numerically the parameter region beside the Toulouse point U = in Section ^. 

We summarize below the advantages to treat the resonance level model rather than the 
Kondo model: First, in the former model, the spin index is dropped so that we can diagonalize 
the system size twice as large as that of the Kondo model. The physical reason of this dropping 
is originated in the foUowings. Because the Kondo problem is essentially of a one-dimensional 
problem. The charge and spin degrees of freedom are separated completely. Therefore, only 
the spin degrees of freedom are subjected to the impurity-spin scattering. That is, the spin 
degrees of freedom transmit through the impurity through the assistance of the scattering, 
whereas the charge degrees of freedom do not transmit at all. That is why we can ignore the 
charge sector. (Therefore, the perfect-transmission conductance is not 2e^//i, which arises in 
conventional models, but /h, because the charge-sector transmission is ignored completely 
here.) Secondly, as is discussed in the previous paragraph, it is very advantageous that we 
can start from the Kondo fixed point, namely, the Toulouse point, which readily contains the 
very essence of the Kondo physics. Hence, we can study the effect of the biquadratic term in a 
systematic manner as the gradual deviation from the fixed point. 



3 Formulation of the nonequilibrium electron transport 



In this section, we discuss how we set up the formahsm to simulate the nonequihbrium electron- 
transport phenomenon. The formahsm is prepared in such a way that it is readily implemented 
in a computer-simulation algorithm. First, we give a conductance formula based on the 'con- 
ductance viewed as transmission' picture proposed by Landauer Q. Extensive use of 
numerical technique enables the evaluation of the transition-matrix elements containing higher 
order contributions of both the hybridization amplitude and the density-density correlation. 
Then, we give the algorithm to calculate the nonequilibrium Green function, which is useful in 
order to investigate the spectral property under nonequilibrium current flow. 

According to the Landauer picture, the electron current is given by the transmission prob- 
ability between two leads; 

/ = eJ2{f{yFk ATkk'Wl- f{vFk'))5{vFk-VFk' + eV) 

kk' ^ 

-(1 - f{vFk))^\Tkk'\^f{vFk')S{vFk - VFk' + eV)^ 

= eN^pl I ds - eV)^\T,_,vJ' (1 - /(e)) - (1 - /(e - eV)) ^\T,.,yJ^ f{e)^ 
r^v/-^ 27r 

= eN^pl / de— |T,_,v-/2,,+ey/2p. (6) 

J-eV/2 " 

We have used the Fermi-distribution function /(e) and the density of states = 1/2, and 
supposed that the system is at zero temperature. Tkk' denotes the transition-matrix element 
from the wave number k of the left lead to fc' of the right lead. The explicit expression of Tkk' 
is considered afterwards. Thereby, we obtain the differential conductance. 



d 9 9/' 27r,„ ,o 



Therefore, the conductance G{eV) is given by the derivative of the integral of the following 
integrand, 

g(e) = eN^p^^^\T,^,v/2.e+eV/2\^- (8) 

(Note that the function g{e) is an even function because of the particle-hole symmetry.) Hence, 
it is expected that the integrand g{e) yields the conductance, 

G{eV) ~ egieV/2), (9) 

unless the integrand depends on eV very much. Fortunately, it is known that in the noninter- 
acting case U = 0, the current is expressed in terms of the following rigorous form |29 



r.'-^ii^r^SMK^ (10) 

" J-eV/2 

(That is, at [/ = 0, the transition matrix Ti;-eV/2,e+eV/2 is related to the local Green function 
at the impurity site (oc Q(uj = e)).) Because the integrand is completely independent on eV, the 
relation ^) holds exactly. We expect that beside U = 0, the relation continues to hold well. In 
fact, as is presented afterwards in Section 4.3, uj{= e)-dependence is dominated and determines 
the essential feature (peak structure) of the integrand, whereas eV just causes sub-dominant 
detailed contributions. As is shown below, a handy formula ( p6| ) for conductance estimate exists 
and is expressed in the similar form as ours (g). In the presence of U, the integrand behaviors 



of ours (^) and the formula become qualitatively different. Therefore, at the present stage, 
it would be well worth considering the conductance in the approximation level of (^) . 

Next, we calculate the transition-matrix element. In order to do that, we must prepare the 
initial and final states, 



and 



J90) 



|f> = 4. l5o) (.? = N/2 + 1, 7V/2 + 2, • • • , N), 



(11) 



(12) 



respectively. We have supposed that N is an even integer. The state \go) represent the situation 
where both leads are half-filled; 



l5o) = 4,4,---4„,,4,<---<,Jo). 



(13) 



In the single-particle transmission process, the energy should be conserved. This conservation 
condition restricts the available values of the voltage drop within the series eV = 2{j — N/2 — 
l)/(iV — 1) {j = N/2 + 1, N/2 + 2, - ■ ■ , N); note that our conduction-electron band is discrete. 
With use of these initial and final states, we obtain the transition probability. 



t=0 



with. 



Ci{t) = (f| 



Ulrj{t, -Oo)|i) 



ic/i^(t,-(X3)|i)r 

In the above, Uir; is the time-evolution operator of the interaction representation; 



f/i^(i2,ti) -Te J'l 



(14) 
(15) 

(16) 
(17) 

(18) 

fc ^ ^ fc ' 

{7i — Tio+Tii.) The normalization factor in eq. (|l5|) is vital, because the infinitesimal damping 
parameter rj violates the unitarity of C/i,,. This fact is known as the Gell-Mann-Low theorem . 
Moreover, it is notable that rj plays a significant role to realize 'nonequilibrium dissipative state' 
breaking the time- reversal symmetry p8[ |. Using the property of the time-evolution operator 
Uiri, we obtain the following expression. 



with. 



and, 



Hiit) 



Wo = ^ ( z;Ffc + ^ ) LlLk + E ( - ^ ) 4^fc- 



eV 



|(f|C/i,(t,-oo)|i)|^ 



dt {i\Ulrf{0O,t)Ulrj{t,~0o)\l) 
d |(f|C/i,(i,-00)|i)|2 



t=0 



dt (i|C/i^(cx),-oo))|i) 



1 



|[/i,(o,-(X3)|i; 



p ^l(f|C/i.(i,-^)|i>l' 



(19) 



Through expanding the operator [/j^ into the Dyson series, one arrives at the following formula. 



2-K 



■\T\i)fS{Ei-EO, 



where the transition matrix T is. 



T : 



1 



|[/i^(0,-oo)|i>| 



Hi 1 



1 



Ho + Hi-Ei- 17] 



-Hi 



(20) 



(21) 



with 7^0 |i> = Ei\i) and, 

C/i.(0,-oo)|i)=(l-^^-^^^^L___^,)|i). (22) 

The above formulae complete our technique to compute the electron conductance; note the 
relation Tq^eV — Readers may have noticed that we are subjected to the Lippmann- 

Schwinger formulation essentially [Q. It is notable that the formula for the transition matrix 
( pl| ) contains infinite-order contributions with respect to Tii, and thus our simulation does not 
resort to any perturbative treatments. 

Lastly, we show the way to calculate the non-equilibrium Green function, 

g{eV,uj) = J dte'"*e(ey,t) (23) 

g{eV,t) - -iQ{t){g,vmS-it),S+}\gev), 

where the nonequilibrium steady state is given by, 

_ [/i^(0,~oo)|go) 
|;7i^(0,-cx))|3o)| 

Because of the same reasoning as the above, we need a normalization factor. The time-evolved 
state is calculated similarly, 

C/i„(0, -^)\go) =(l- ^ \ -Hi] |.go). (25) 

This gives an explicit expression for the state \gev)- What is still left is to calculate the Fourier 
transform of the time correlation function of eq. (p3|). At first glance, one might think that it is 
impossible, because all the eigenstates are needed to evolve the time correlation. Gagliano and 
Balseiro, however, invented a way to express the Green function (j2^) in a compact continued- 
fraction form with the Lanczos tri-diagonal elements p9[ . Through utilizing their technique, 
the scheme to calculate nonequilibrium Green function is completed. We emphasize that in our 
scheme, we do not resort to any perturbative treatments. 

The Green function is useful to estimate approximate value of electron current. So far, the 
following formula has been used to obtain the nonequilibrium current; 



e 



V/2 pLpJ?, _-| 

,^_^^^^d.j^.-lmg(ey,.), (26) 

with the Abrikosov-Suhl resonance frequency for each lead namely, 

rL.i^ = 2^p,(K/\/2)^ (27) 



As is mentioned in Introduction, this formula ( p6D insists that the spectral density at the 
impurity site (oc ImQ) contributes to the current flow. Note that as in our above-mentioned 
formula (^, the current is expressed in terms of an integral form. An so, because of the same 
reasoning mentioned before, the integrand of (p6| ) does hardly depend on eV for J7 w 0, yielding 
a compact expression for the differential conductance; 



g2 pLpi? 



In the next section, we report that the conductance estimate with our formula (^) differs 
qualitatively from that of the above handy relation. 



4 Numerical results and discussions 



So far, we have prepared prescriptions for simulating the nonequilibrium electron transport 
numerically. In this section, we carry out numerical simulations based on the prescriptions. 



4.1 Preliminaries of our numerical simulation 



Here, we summarize details of our numerical computation. We simulated the system consisting 
of twenty conduction-electron levels. (The form of the Hamiltonian is explained in Section |2f 
see eq. (|l|).) That is, there are ten levels for each lead; N = 10. (Owing to the presence of the 
density-density-correlation term in the Hamiltonian (|l|), the Hamiltonian- matrix elements are 
not so sparse as in conventional diagonalization simulations. That costs considerable compu- 
tation time, even though the dimensionality of the Hilbert space (~ 2^^) is not so particularly 
overwhelming . ) 

The most important computation stage is that to evaluate the resolvent which is appearing 
in eqs. (|2|) and (|2^). The resolvent is computed with the conjugate-gradient algorithm. The 
convergence of the conjugate-gradient iteration is not very stable. This instability becomes 
more serious as the coupling to the leads is strengthened. This instability has been encountered 
so far in general in calculating resolvent of many-body problem numerically. Here, we have 
used the following new trick to overcome this difficulty We substituted the parameter rj in eqs. 



(y) and (122|) with, 

rj = r,(") + r;(2) ^((^Ffc)^^^^ + {vFk)^RlRk). (29) 

k 

Now, 77 is not a c-number constant, but is an operator. Because the parameter determines the 
energy resolution, it should be of the order of the conduction-band discrete spacing. Hence, we 
set T]^^^ =0.1. The first term alone results in a desperate instability of numerical procedure. 
Hence, we have added the second new term. The meaning of this term may be transparent: The 
term enforces the life time of each conduction-band level to vary as the square of the excitation 
energy measured from the Fermi level. Such dependence might be reasonable, because the 
Landau theory for the interacting electron gas concludes the same dependence of the life time. 
Here, we set 77^2) = 0.06 for Vh = 0.3. 

We found that this choice of parameters is the best one: Reader may feel that our N = 
10-level approximation of the conduction band might be serious, because the Fermi-surface 
singularity is important to realize the Kondo effect. Yet, owing to the resonance-level-model 
mapping, we are actually in the stronf^-Kondo-coupling regime (Toulouse point), where the 
Kondo temperature is given by ~ F'^ -I- F^ = 271 p^V^; see eq. (p7|). As a matter of fact, our 
choice of Vh = 0.3 gives the Kondo temperature ^ 0.28 which is larger than the conduction- 
band level spacing 2/N = 0.2. On the contrary, the energy-resolution broadening term 77 is not 
negligible compared with the Kondo temperature so that we are suffered from a broadening of 
the zero-bias-anomaly peak. We found that for Vh > 0.3, the numerical instability becomes 
desperately serious, and exceedingly large smearing factor rj is needed. Hence, we concentrated 
ourselves on the above-mentioned parameter point with Vh — 0.3, Ty'^"-' = 0.1 and 77^^^ = 0.06, 
where we can achieve the Kondo effect and barely manage the numerical instability with a 
modest energy-resolution broadening. 



4.2 Differential conductance 

In this subsection, we present our numerical result of the differential conductance. The calcu- 
lation is based on the formula (||). We discuss the result with an emphasis on the effect of U, 
namely, the deviation from the Toulouse point. 

In Fig. 1^, we plotted the differential conductance against the voltage drop eV for Vh ~ 0.3 
and various U. We see that the differential conductance is maximal at the zero bias eV = 
irrespective of U. This peak structure is known as the 'zero-bias anomaly' The anomaly 
has been observed experimentally, and is very significant to confirm that the electron actually 
tunnels between two leads with a certain quantum-mechanical tunneling amplitude. 

First, we discuss the data for i7 = (the Toulouse point). At this quadratic case [/ = 0, an 
analytical prediction for the differential conductance is available; see Section 0. According to 



that, the differential conductance is given by the formula (Pq ) with the Green function, 



g{u;) 



1 



(30) 



k w+S-VFk-eV/2 



i^+S-VFk+eV/2 



The result of this formula with 5 = 0.05 is also plotted as a dashed line in Fig. We see 
that our numerical result succeeds in reproducing the zero-bias conductance ~ e^//i, which 
is accordant with the Landauer theory, and the zero-bias peak as well. On the contrary, our 
zero-bias-peak width is broadened substantially compared with the analytical prediction curve. 
That is due to the energy-resolution-smearing factor rj (p^), which cannot be omitted in order to 
stabilize the conjugate-gradient procedure. (This factor 77 broadens the energy resolution, and 
loosens the energy-conservation constraint. Thereby, this parameter enhances the transmission 
probability, resulting in an enhanced estimate of conductance.) As is explicated in Section 4.1, 
we had implemented a new trick that the factor 77 (life time) is not a constant, but is gained 
near the band edges. Therefore, owing to 77, our conductance data suffer a correction and are 
enhanced especially for large potential drop eV ^ 1. Hence, we rather focus our attention 
on the qualitative variation of the conductance curve G{eV) due to the introduction of the 
many-body correlation U. It should be stressed that, as is explained in Section ^ besides this 
broadening, our simulation is free from any approximations, and actually, takes into account 
full many-body correlation processes. 

Second, keeping that in mind, let us turn to the case of J7 < 0. Note that the case U < 
corresponds to the situation where the Kondo coupling is reduced from the value of the Toulouse 
point; refer to the mapping relations (^). In Fig. |^, we see that the zero-bias conductance is 
suppressed significantly by the reduction of the Kondo coupling. It is very surprising, because 
conventional picture based on the one-particle description, namely, eqs. (^8|) and (|2^), gives 
the conductance e'^/h irrespective of the hybridization amplitude as far as the left-right hy- 
bridization couplings are symmetric. Therefore, we notice that such reduction is of a nontrivial 
many-body effect, and in principle, that should be managed by ab initio treatments such as 
ours. As a matter of fact, a conductance estimate based on the convenient formula ( p8|) gives 
qualitatively opposite behavior that the conductance is enhanced with U < 0; this discrepancy 
is reported in the next subsection. 

Furthermore, it seems that the width of the zero-bias-anomaly peak is not changed very 
much by [/ < 0. It may be interesting to compare our findings with those of the lowest-order- 
perturbation theory with respect to the Toulouse point by Majumdar, Schiller and Hershfield 
[p^ . They studied the effect of certain three types of perturbations to the differential con- 
ductance. These perturbations are thought to drive the integrable Toulouse Hamiltonian in 
the direction of the weak Kondo coupling; namely, more realistic coupling. (Because they are 
employing a field-theoretical description together with a very sophisticated canonical transfor- 
mation, the origin of these perturbations in the language of the original Kondo Hamiltonian 
is rather unclear, and is to be clarified.) They found that some two perturbations keep the 
height fixed, but sharpen the width, whereas the other perturbation suppresses the height with 
the width unchanged. The true effect is speculated to be a certain mixture of those effects. 
Our results suggest that the reduction of the Kondo coupling reduces the height, but does not 
change the width very much; that is, in the field-theoretical language |p^, the latter type of 
perturbation might be realized actually in the physics of [/ < 0. 

Next, we discuss the case U > 0. This repulsive density-density correlation is interpreted 
as the very strong Kondo coupling, which is even exceeding the Toulouse-coupling strength. 
Hence, this case might not be so relevant in the language of the Kondo physics, and has not yet 
been studied very well |Q. Nevertheless, apart from the context of the Kondo problem, this 
situation is interesting by its own right. In Fig. |l], We observe that as the repulsive correlation 
is increased, no particular change is observed for the differential conductance especially for 
eV 0. That is, despite of the repulsion, inherent quantum-mechanical resonance restores the 
hybridization between the impurity and the conduction electron. This situation will be more 
clarified in the next subsection through referring to the spectral-function data. 



4.3 Nonequilibrium Green function 



In this subsection, we present the numerical resuh of the nonequihbrium Green function. The 
calculation is based on the formula ( p3|) . We concentrate on the imaginary part of the Green 
function, which yields the local density of states at the impurity site; 

p{uj) = -llmg{eV,uj + iS), (31) 

TT 

with 5 = 0.1. Once the spectral function is obtained, one can adopt the handy formula ( p^ ) 
to estimate the conductance. We show that the estimates (behaviors) differ qualitatively from 
our first-principle results shown in Section 4.2. 

In Figs, we plotted the density of states for Vh = 0.3, and U — 0, 0.2 and —0.2, 

respectively. The local density of states shows a peak structure. That is, the impurity level 
is smeared out by the mixing (resonance) to the conduction-electron band. The peak width 
corresponds to the inverse of the life time of the impurity state. 

First, we discuss the case of i7 = (Fig. ||). We see that the local density of states, 
namely, the resonance peak, is suppressed gradually as the bias voltage eV is increased. That 
is understood as follows: For eV ~ a;c(= 1), the Fermi-energy position approaches the band 
edge. In that case, the resonance is not fully formed because of the disturbance of the band 
edge. That causes the reduction of the resonance peak for eV ~ 1. 

Second, let us turn to the cases with the many-body correlation U ^0. In those cases, the 
Green function becomes far more nontrivial, and contains significant informations. For U = 0.2 
(Fig. ^), the density of states forms broader resonance than that of U = 0. In particular, the 
sub-peaks a,t uj — zLeV/2 are prominent; note that the frequencies uj — zLeV/2 are the Fermi 
levels of the leads. The sub-peak height stays aloft even for eV > 0. Thereby, the sub-peaks 
constitute shoulders of the main peak, resulting in considerable broadening of the main peak. 
This feature has been captured by previous studies |]l^ [l^, ^ , and is speculated to be the very 
essence of the nonequilibrium Kondo physics. We emphasize that the present simulation does 
actually capture this characteristic without resorting to any approximations. This prominent 
resonance may be the precursor of the unexpected large conductance, which we reported in the 
previous subsection. On the other hand, for U = —0.2 (Fig. ||), the density of states forms a 
very narrow peak, and the shoulders at w = ±6^/2 are not grown very much. That is, in this 
weak Kondo couphng case, the Kondo resonance is suppressed to a considerable extent. This 
is consistent with the observation of the previous subsection that the (zero-bias) conductance 
is suppressed significantly hy U < 0. 

Finally, we show the result of the differential conductance based on the approximative 
formula (|2^). The application of the formula is now possible, because we have nonequilibrium- 
Green-function data. The result is plotted in Fig. ^ Although this approximative formula 
reproduces the zero-bias anomaly, the dependence on U differs significantly from that of our 
first-principle result. For instance, in Fig. |^, the conductance is increased by i7 < 0, whereas 
it should be suppressed according to our simulation. (Our conclusion may also be supported 
by a recent report, where, however, the authors studied the tunneling amplitude from one edge 
to the bulk of the Hubbard chain 0].) It is suspected that the approximate formula does not 
fully include the many-body correlation effect. One cannot expect, as a matter of fact, that 
the nonequilibrium transport coefficient is given by one-particle Green function alone, because 
the Kubo theory tells that the 77.on-linear-response function should be expressed in terms of 
many-point Green functions. 

5 Summary 

We have investigated numerically the nonequilibrium electron transport through the Kondo 
impurity. In this problem, nonequilibrium-driving force and biquadratic correlation coexist, 
and both are playing a crucial role. For the first time, we succeeded in simulating the situation 
without resorting to any perturbative treatments. Moreover, because we utilized the mapping to 



the resonance level model, we could clarify the effect of the biquadratic correlation as a gradual 
deviation from the Toulouse point. We calculated the differential conductance with the formula 
and the nonequilibrium Green function with eq. (p^). Our results tell that the effect of 
the biquadratic interaction lies out of the scope of the one-particle description: The attractive 
density-density coupling, which corresponds to a reduced Kondo coupling from the Toulouse 
point, suppresses the zero-bias conductance below the Landauer value. The width of the zero- 
bias anomaly, on the other hand, is not changed very much. The above is to be contrasted 
with the field-theoretical lowest-order-perturbation theory with respect to the Toulouse point 
. The repulsive coupling, on the contrary, does not influence the conductance. The above 
features are supported by our spectral function data based on eq. (|3l|). 

We stress that as we have demonstrated, the 'steady nonequilibrium' transport can be 
simulated without carrying out any time-evolution simulations. The formalism and the imple- 
mentation to the computer algorithm, which are demonstrated here, are readily applicable to 
other wide class of steady nonequilibrium transport with many-body interaction. As is clar- 
ified here, in such situation as well, the many-body correlation would cause new unexpected 
behavior. This would be remained in future. 
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Figure 1: Differential conductance is plotted for Vh = 0.3 and various U . The dashed curve 
shows the result with the analytic formula ( ^ ) that is valid at [/ = 0. 




Figure 2: Local density of states at the impurity site ( |3l| ) for Vh — 0.3 and U = 0. 




Figure 3: Local density of states at the impurity site (^) for V|i = 0.3 and U = 0.2. 




Figure 4: Local density of states at the impurity site (|3l| ) for V\i — 0.3 and U = 



-0.2. 




Figure 5: Differential conductance evaluated by means of the approximate formula (28) with 
our spectral-function data shown in Fig. This is to be compared with the first-principle 

result shown in Fig. |^. 



